function solution = wrap_PATH( obj )
    A = obj.dynamics.A_MCP;
    b = obj.dynamics.b_MCP;
    z0 = obj.dynamics.z0_MCP; 
     
    nb = length(obj.dynamics.Vel) / 6;
    nc = length(obj.contacts.PSI);
    nd = obj.dynamics.nd;
    U = obj.dynamics.U;
    PSI = obj.contacts.PSI';
    
    problem_size = length(b);
    big = 10^20;
    u = big*ones(problem_size, 1);
    l = [-big * ones(6*nb, 1);
          zeros((2+nd)*nc, 1)];
    obj.solution.big = big;
    obj.solution.l = l;
    obj.solution.u = u;
    tic;
    solution = pathlcp(obj);
    solution.solve_time = toc; 
    obj.z = solution.NUnew;
%     error = z' * (A * z  + b);
%     solution.z = z;
%     solution.iter = 0;
%     solution.err = error / nc;
%     
%     % the detailed error information
%     Gn = obj.dynamics.Gn;      
%     Gf = obj.dynamics.Gf;
%     U = obj.dynamics.U;
%     
%     NU = z(1:6*nb, 1);
%    
%     rel_vel_n = Gn' * NU;     % nc X 1
%     rel_vel_t = Gf' * NU;     % (nc*nd) X 1   
%     % normal error
%     norm_err = pn' * rel_vel_n;
%     % resultant frictional force and velocity on the tangential plane
%     resultant_f  = zeros(nc, 2);
%     resultant_v  = zeros(nc, 2);
%     theta = 2*pi/nd;
%     % get the friction and velocity in the tangential plane as resultant force and velocity. 
%     slide_err_angle = 0;
%     slide_err_mag = 0;
%     for i = 1 : nc
%         for  j = 1 : nd
%             resultant_f(i, 1) = resultant_f(i, 1) + pf((i-1)*nd+j, 1) * cos((j-1)*theta);
%             resultant_f(i, 2) = resultant_f(i, 2) + pf((i-1)*nd+j, 1) * sin((j-1)*theta);
%             resultant_v(i, 1) = resultant_v(i, 1) + rel_vel_t((i-1)*nd+j, 1) * cos((j-1)*theta);
%             resultant_v(i, 2) = resultant_v(i, 2) + rel_vel_t((i-1)*nd+j, 1) * sin((j-1)*theta);
%         end八零
%     end
%     for i = 1 : nc
%         vel_mag = sqrt(resultant_v(i, 1)^2 + resultant_v(i, 2)^2);
%         fric_mag = sqrt(resultant_f(i, 1)^2 + resultant_f(i, 2)^2);
%         if vel_mag > 1e-6 && fric_mag > 1e-6
%           % sliding case
%           vel_dir = resultant_v(i, :) / vel_mag;
%           fric_dir = resultant_f(i, :) / fric_mag;
%           slide_err_angle = (1/pi)*(acos(vel_dir*fric_dir') - pi);
%           slide_err_mag   = norm(U(i, i) * pn(i, 1) - fric_mag );
%         else
%           % sticking case
%           stick_err = min(U(i, i) * pn(i, 1) - fric_mag, 0); 
%         end
%     end
%     solution.norm_err = norm_err /nc;
%     solution.fric_err = (slide_err_angle + slide_err_mag + stick_err) / nc;
%     solution.slide_err_angle = slide_err_angle;
%     solution.slide_err_mag = slide_err_mag;
%     solution.stick_err = stick_err;
end
